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We present a new exact algorithm for estimating all elements of the quark propagator. The advantage of the 
method is that the exact all-to-all propagator is reproduced in a large but finite number of inversions. The efficacy 
of the algorithm is tested in Monte Carlo simulations of Wilson quarks in quenched QCD. Applications that are 
difficult to probe with point propagators are discussed. 



1. INTRODUCTION 

The discovery of exactly chirally symmetric 
fermions on the lattice has triggered intensive re- 
search in the development of algorithms to simu- 
late dynamical Ginsparg- Wilson fermions on the 
lattice. In light of recent work in this area, we 
anticipate that there will be a limited number of 
expensive, full QCD configurations in the near 
future. One would like to extract all the informa- 
tion that one can from these lattices without be- 
ing restricted by point propagators, which would 
appear highly wasteful considering the cost to 
generate the configurations. Point propagators 
do not require massive computing power but re- 
strict the physics one has access to, mainly the 
flavour non-singlet spectrum. They also restrict 
the interpolating operator basis used to produce 
early plateaux in effective masses, for instance, 
since a new inversion must be performed for every 
operator that is not restricted to a single lattice 
point. Variational methods would be much more 
powerful with the use of all-to-all propagators. 

All-to-all propagators |H2l3l4IE>lrj] provide a so- 
lution to these problems, but are usually too ex- 
pensive to compute exactly as this requires an un- 
realistic number of quark inversions. Stochastic 
estimates tend to be very noisy and variance re- 
duction techniques are crucial in order to separate 
the signal from the noise. In this paper we pro- 
pose an exact algorithm to compute the all-to-all 
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propagator utilizing the idea of low-mode domi- 
nance corrected by a stochastic estimator which 
yields the exact all-to-all propagator in a finite 
number of quark inversions. 

2. THEORY AND NOTATION 

2.1. Spectral Decomposition 

Theoretical arguments, backed up by numer- 
ical evidence |7I8| . indicate that the low lying 
eigenmodes of the Hermitian Dirac operator cap- 
tures much of the important infrared physics in 
hadronic interactions. It would be desirable to 
take advantage of this fact and solve exactly for 
as many of the eigenmodes as possible to estimate 
the all-to-all quark propagator. 

To this end, we define the Hermitian Dirac op- 
erator Q = 75 M where M is the usual Dirac oper- 
ator. The truncated representation of the quark 
propagator is then given by, 

Oo^r^Ki/o)®^^^, (1) 
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where Qv^ 1 ' = \iV^\ The truncation at an ar- 
bitrary number, N ev , of eigenvectors, however, 
leaves the theory non-unitary, making it manda- 
tory to correct it. 

2.2. Dilution (Noisy Estimators) 

The standard method of estimating the all- 
to-all quark propagator is to sample the vector 
space stochastically. One generates an ensemble 
of random noise vectors, {wii, - '' ,V[N r ]}' w hh 
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the property 

((ri(x)®ri(yy)) cx 8 x>y , 



(2) 



where ((• • •)) denotes average over noise samples. 
The solution for each vector is obtained in the 
usual way, 



i>[r] ( x ) = Q 1? ?[r] (y) ■ 



(3) 



The all-to-all quark propagator is estimated as 

Q- 1 = U®^)) 



(4) 



This method is noisy because it relies on delicate 
cancellations in the 0(1) noise over many samples 
to find the signal, which falls off exponentially 
with the separation. We propose to remove the 
0(1) random noise by "diluting" the noise vector 
according to some dilution scheme resulting in an 
exponential gain in the variance. A particularly 
important example of dilution for measuring tem- 
poral correlations in hadronic quantities is "time 
dilution" where the noise vector is broken up into 
pieces which only have support on a single times- 
lice each, 
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where rf^> (x, t) — unless t = j. 

Each diluted source is inverted resulting in N t 
pairs of vectors, {ijjW (x, t),r]^' (x, t)}, which then 
gives an unbiased estimator of (ipip) with a single 
noise source, 
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We show the effect of time dilution on a pseu- 
doscalar propagator on a 12 3 x 24 lattice in Fig.^ 
The triangles are the average of 24 noise sources 
without any dilution and the circles are from a 
single noise source which has been time-diluted. 
This particular scenario is analogous to the "wall 
source on every timeslice" method used by the 
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Figure 1. The psuedoscalar propagator computed 
with and without time dilution. 



authors of Ref. pQ to estimate the disconnected 
diagrams appearing in hadronic scattering length 
calculations. Our method is, however, more gen- 
eral and can be extended to the spin, color and 
space components of the source vector. The 
"homeopathic" limit of this dilution procedure re- 
sults in the exact all-to-all propagator in a finite 
number of steps (Fig. 0). This limit cannot be 
reached in practice on realistic lattices, but the 
path of dilution may be optimized so that the 
noise from the gauge fields dominate the errors 
in the hadronic quantities of interest with only a 
small, managable number of fermion matrix in- 
versions. 

2.3. 'Hybrid' Method 

The evidence that much of the hadronic physics 
we would like to extract from the lattice is in the 
low-lying eigenmodes of the (Hermitian) Dirac 
operator suggests that one try to calculate as 
many as possible of the low modes exactly and 
correct for the truncation with the noisy method. 
This gives rise to two concerns: firstly, the cor- 
rection should leave the exactly solved low-lying 
modes intact; and secondly, it should not intro- 
duce large uncertanties in the process. We pro- 
pose that the stochastic method with noise dilu- 
tion is a natural way to accomodate both of those 
concerns. 

First, we note that the exact N ev low eigen- 
modes obtained separately naturally divide Q 
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Figure 2. A cartoon of possible deviations of 
the stochastic estimates of the exact solution (at 
Nd = Nraax) for different dilution paths. Simply 
adding noise vectors will give a \/\f~N behaviour. 



into two subspaces, Q = Qo + Qi, defined by 



Qo = ]TW 4) <W 4)t , 
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Similarly, the quark propagator is broken up into 
two pieces, Q _1 = Qo + Qi, where Qo is given by 
Eq. Q and Qi is the remaining unknown piece. 
We correct for the truncation and estimate Qi 
using the stochastic method, Qi = ({ipi&rf}) with 
N r noise vectors, {?7[i], • • • , V[N r ]}- The solutions 
are given by 

VV] = {Q^Vl) V[r] = (VlT] [r] ) , (9) 

where V\ is the projection operator 

N ev 

V 1 = l-V = l-J2 vU) ® w0)t ■ ( 10 ) 



Note that Eq. © follows fr om the identity 
Q _1, Pi = Qi. As was mentioned earlier, the idea 



of dilution will be applied to the stochastic es- 
timation of Qi. Each random noise vector, n r ], 
that is generated will be diluted and orthogonal- 
ized (wrt «W) so that it can be used to obtain 
ip\ r ]. In other words, we now have the following 
set of noise vectors: 
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where the upper indices denote the dilution and 
the lower indices label the different noise sam- 
ples. We note that the noise vectors are mutually 
orthogonal to each other due to the dilution be- 
fore an average over different random vectors are 
taken, i.e., 



(11) 



This results in smaller variance than the standard 
method which mixes noise, as Eq. shows. 

There is a natural way to combine the two 
methods to estimate the all-to-all quark propa- 
gator. The similarity in the structure of Eq. 
and Eq. J3J suggests that one construct the fol- 
lowing "hybrid list" for the source and solution 
vectors: 
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where the indices run over Nhl — N ev + Nd el- 
ements. The master formula for the unbiased, 
variance reduced estimate of the all-to-all quark 
propagator is then given by 



Nhl 



M~ 



= (x, x A ) <g> w {l) (y, y 4 ) f 75 • (14) 



Using the pion as an example, we can demon- 
strate how unitarity is recovered from the trun- 
cated propagator. In Fig. |3 we show the effective 
mass from the truncated propagator and from 
the hybrid method with a time, spin, color and 
space (even-odd) diluted noise vector. The trun- 
cated propagator, which approaches the asymp- 
totic value from below, is corrected by the addi- 
tion of the diluted noisy propagator. 
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Figure 3. The effective mass from 100 eigenvec- 
tors and from the hybrid method with 100 eigen- 
vectors and a time-diluted noise vector. The lat- 
tice size is 6 3 x 12. 



3. QCD IMPLEMENTATION 

All-to-all propagators not only expand the 
range of applications accessible in lattice QCD 
but also considerably simplify the construction of 
hadronic interpolating operators. A pseudoscalar 
correlator is constructed in the following way with 
traditional, point-to-all quark propagators, 

C = M- x {S,t; 0,0)^ M^&tAO)^. (15) 

The construction of a better operator is awkward 
at best since the operator at timeslice t is con- 
structed from various components of the quark 
propagator connecting the source and sink. All- 
to-all propagators eliminate this complication as 
the operator at timeslice t is constructed from 
vectors on that timeslice, 



(16) 



where the extra factor of 75 comes from the use 
of the Hcrmitian Dirac matrix, 75M. For compli- 
cated operators such as those used to project out 



hybrid/exotic states, this is a much needed sim- 
plification. For example, an interpolating opera- 
tor for the exotic hybrid 1 h can be constructed 
from combinations of gluonic paths projecting out 
the relevant quantum numbers [21 ■ One term in 
such a sum may be 

w^(x^U z {x)U y {x+z)Ut(x+ey)u [j] {x+e y ) , (17) 

where all the variables are located on a single 
timeslice. A standard P-wave state may for ex- 
ample be constructed as follows, 



0$ j} (x, t) = W W (x, i)t(D fc utf>)(af, t) , 



(18) 



where Dk is the covariant derivative. 

Using all-to-all propagators, correlation func- 
tions are also constructed in an intuitive man- 
ner. Hadronic correlation functions are obtained 
by correlating interpolating operators sitting at 
different timeslices, e.g. 
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for isovector two-point correlators (propagators). 

Noting that this is simply a contraction of the 
source and solution vectors with an outer sum 
over the hybrid list indices, this gives even greater 
scope for simplification. Once the machinery of 
the hybrid list summation is set up in the pro- 
gram, it becomes a black box to the end user who 
simply supplies the subroutines to create the re- 
quired operators from the quark, antiquark and 
gluon fields on a timeslice. 

For isoscalar mesons, the disconnected part of 
the propagator must be included, yielding the fol- 
lowing contraction, 

-{^WSr^V)} { w «(io) t 7 5 r t w g(to)} • 

4. TESTS 

We use the Wilson action to illustrate the ef- 
fectiveness of the variance reduction although it 
is expected to work even better for a chiral action 
and light quarks. The simulation parameters are 
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Table 1 

Fermion run parameters, 
are taken from Ref. 



The values for m. 
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Figure 4. Equal cost test of the pion effective 
mass at timeslice 6. 



shown in Table ^ In order to demonstrate the 
efficacy of the noise dilution method, we perform 
an "equal cost" comparison test. Equal cost here 
shall mean the same number of inversions per- 
formed, disregarding the sparseness of the source 
vector. For example, an even-odd dilution in 
space results in twice as many inversions as the 
one without, hence two independent sets of time 
diluted noise vectors are generated to compare 
the noise in the effective mass of the psuedoscalar. 
As we have seen that no signal is obtained oth- 
erwise, all noise vectors have been time diluted. 
The results are shown in Fig.0J The noise diluted 
errors consistently show smaller errors than the 
standard stochastic method. At this point, one 
may note that the isovector correlation function 
(when both quarks are saved as all-to-all) involves 
a large number (Nhl x Nhl) of contractions. We 
can take advantage of having saved different ran- 
dom source/solution samples by reusing them in 



the contraction. In other words, one can gener- 
ate Nr samples of noise vectors, r)[ r ], and save the 
corresponding solutions, ^>r r i to disk and perform 
the contraction, 



C 



u m(*W s i. u m(*o) ■ 



(20) 



yielding ~ N R samples of the pion correlation 
function. The errors correspondingly decrease 
faster than the naive l/V Nr, although the mea- 
surements are correlated. We have seen in our 
preliminary tests that the gain in error reduction 
is comparable to some dilution choices. It is clear 
that if one can afford to save the noise/solution 
vectors onto disk, then this is a straightforward 
method of variance reduction. 

We have performed some spectroscopy calcu- 
lations using the idea of dilution. Standard me- 
son spectroscopy including pions, rhos, heavy ex- 
otic hybrids (with two different operators) and 
standard P-wave states were studied. For the P- 
wave states in particular, it is difficult to get a 
good overlap using point propagators. The ef- 
fective masses for all the particles are shown in 
Fig. O and Fig. Although in this test study 
we only used 10 quenched gauge configurations, 
we have a surprisingly good signal for P-waves as 
well as the extremely heavy hybrid states. Fur- 
thermore, decay constants can be extracted from 
point-smeared and smeared-smeared correlators 
if all-to-all propagators are used ^U] ■ 

5. SUMMARY 

We have presented a new algorithm to estimate 
the all-to-all propagator. All-to-all propagators 
make it possible to make use of all the available 
information in a gauge configuration, which con- 
sidering the cost to generate full QCD configura- 
tions may be of crucial importance. They are also 
a necessary ingredient in flavour singlet physics. 

All-to-all propagators have a further advantage 
over point propagators in that operator construc- 
tion is considerably simplified: The operators are 
constructed in a natural way from local fields, and 
extended operators used in variational methods 
may be employed at no additional cost. 

We have presented evidence that diluting 
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Figure 5. The pion and rho meson effective 
masses from 10 configurations on a (3 — 5.7, 
1 2 3 x 24 lattice. 



Figure 6. Preliminary results for a variety of P- 
wave states. The staple hybrid operators actually 
use a length-2 staple. 



stochastic estimators in time, colour, spin or 
other variables gives less variance than traditional 
noisy estimators. This is not unexpected, since 
dilution will yield the exact all-to-all propagator 
in a finite number of iterations. More work is 
needed to determine the optimal dilution path for 
different observables. 

The hybrid method allows one to extract the 
important physics from low-lying eigenmodes and 
combine this with a noisy correction in a natural 
way. This may be implemented in such a way 
that the user can be blind to the details of the 
dilution and eigenvector list. 

Preliminary tests of the algorithm have been 
presented with small lattices and a small num- 
ber of configurations. We find that the algo- 
rithm works as expected for standard light-light, 
static-light and exotic meson spectroscopy. Fur- 
ther work will include higher statstics on lighter 
masses as well as tests on disconnected diagrams. 
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